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ABSTRACT 


Sauer  recently  proposed  a  new  approach  to  hyperbolic  problems  Involving 
more  than  two  Independent  variables.  In  the  case  of  plane,  unsteady  problems 
depending  on  time  t  and  Cartesian  coordinates  x,y  he  finds  the  characteristic 
lines  of  the  governing  equations  In  the  x.t  plane  and  develops  a  finite 
difference  scheme  based  on  these  lines  and  lines  parallel  to  the  y  axis. 

In  this  way  a  two  dimensional  unsteady  problem  Is  solved  as  a  sequence  of 
one  dimensional  problems  In  planes  y  ■  constant.  This  scheme  Is  much  simpler 
than  earlier  schemes  based  on  bicharacteristics .  The  present  report  describes 
the  scheme  and  discusses  its  application  to  underwater  shock  wave  problems. 


1.  INTRODUCTION 


Many  important  problems  connected  with  the  theory  of  underwater 
explosions  depend  on  two  space  variables  and  time.  In  their  non-linear 
formulation  these  can  be  solved  numerically  by  the  method  of  characteristics. 

When  more  than  two  independent  variables  arise,  the  method  of  characteristics 
can  be  presented  in  a  variety  of  forms  and  it  appears  now  that  some  of  the  early 
numerical  schemes  proposed  for  the  method  were  unnecessarily  complicated.  The 
present  report  describes  some  schemes  which  have  recently  been  developed  in 
Germany,  with  particular  emphasis  on  Sauer's  Near  Characteristics  Method 
(ref.  1),  which  seems  to  be  well  suited  for  application  to  underwater  problems. 
After  considering  the  method  in  its  general  aspects,  it  is  worked  out  in  detail 
for  application  to  two  problems:  (1)  the  propagation  of  a  spherical  explosion 
in  an  ocean  with  a  variable  density  gradient;  (2)  the  reflection  and  refraction 
of  a  spherical  explosive  wave  at  the  ocean  surface. 

Methods  of  Characteristics  developed  earlier  for  two  dimensional  unsteady 
flow  problem  were  based  wholly  on  the  use  of  bicharacteristics.  Those  proposed 
by  Thornhill  (ref.  2)  and  Butler  (ref.  3)  were  based  on  the  following  approach. 
Starting  with  the  Eulerlan  equations  for  unsteady  flow  the  equations  defining 
bicharacteristics  are  determined,  together  with  the  compatibility  conditions 
satisfied  along  them.  It  is  then  assumed  that  initial  data  are  prescribed  on 
some  space-like  initial  surface.  To  construct  conditions  at  a  new  point, 
bicharacteristic  lines  are  drawn  backwards  through  the  point  to  Inteaect  the 
initial  surface.  Along  each  of  these  bicharacteristic  lines  the  compatibility 
conditions  are  written  as  difference  relations.  These  are  simultaneous  equations 
to  determine  the  values  of  the  unknowns  at  the  new  point.  As  many  bicharacteristic 
relations  are  used  as  these  are  unknown  dependent  variables.  Butler  has  applied 
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his  method  to  calculate  the  hypersonic  flow  past  a  pointed  unsymmetrical  body 
(using  the  analogy  between  hypersonic  small  disturbance  theory  and  unsteady 
flow  theory).  The  calculation  was  a  notable  achievement,  although  the  method 
is  more  complicated  than  those  developed  later.  The  chief  disadvantage  of  the 
method  is  that  the  directions  along  which  difference  relations  are  established 
follow  no  particular  pattern.  Further,  it  is  questionable  whether  bicharacter¬ 
istic  relations  on  Mach  cones  can  be  used  more  than  twice  between  an  initial 
surface  and  a  new  point. 

Methods  were  developed  for  steady  flow  problems  by  Coburn  and  Dolph 
(ref.  4)  and  Holt  (ref.  5)  which  only  employ  two  bicharacteristic  relations 
at  each  point.  These,  together  with  a  third  non-characteristic  relation  and 
the  streamline  direction,  are  related  in  a  definite  way  to  prescribed  conditions 
on  an  initial  surface. 

The  methods  developed  recently  in  Germany  also  use  only  two  bicharacter- 
istic  directions  (or  in  the  case  of  Sauer's  method  two  general  sections  of 
the  local  Mach  conoid)  but  the  choice  of  these  is  related  to  the  local  flow 
geometry  rather  than  to  initial  conditions.  The  first  of  these  methods  was 
proposed  by  Bruhn  and  Haack  (ref.  6)  and  later  modified  by  Schaetz  (ref.  7). 

The  basic  approach  is  to  replace  the  derivatives  occurring  in  the  equations 
of  motion  with  respect  to  the  original  independent  variables  (for  example  two 
Cartesian  coordinates  and  time)  by  directional  derivatives.  These  are  chosen 
so  that,  in  the  transformed  equations,  the  number  of  directions  is  reduced  to 
a  minimum.  When  these  transformed  equations  are  replaced  by  difference 
equations,  the  latter  can  be  solved  with  the  minimum  of  algebraic  operations 
and  hence  with  the  least  expense  of  programming  time. 

Xn  the  method  developed  by  Schaets  for  unsteady  flow  the  equations  of 
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motion  are  transformed  to  coordinates  based  on  the  following  local  directions: 

(1)  the  direction  of  the  resultant  velocity. 

(2)  the  direction  of  the  particle  path. 

(3)  bicharacteristic  directions  perpendicular  to  the  resultant 
velocity  direction. 

Schaetz  calls  his  method  optimal  since  it  employs  the  minimum  number  of 
directions.  It  is  a  much  simpler  scheme  than  those  based  entirely  on 
bicharacteristics  -  of  the  four  directions  used  three  are  in  the  same  plane. 
However,  these  directions  do  change  with  the  flow  direction.  Sauer's  Near 
Characteristics  Method  (ref.  1)  avoids  this  drawback.  It  is  also  an  optimal 
method  (using  only  four  directions)  but  uses  a  fixed  coordinate  direction 
as  one  of  these.  The  plane  containing  the  other  three  directions  is  always 
in  a  fixed  direction  (normal  to  the  fixed  coordinate  direction),  two  directions 
are  sections  of  the  local  Mach  cone  on  this  plane,  the  third  is  the  projection 
of  the  velocity  vector  on  the  plane. 

Using  Sauer's  approach  the  method  of  characteristics  for  two 
dimensional  unsteady  problems  reduces  to  applying  the  one -dimensional  method 
in  a  sequence  of  coordinate  planes.  This  makes  it  easier  to  control  than  methods 
based  on  many  spatial  directions  and  easier  to  use  for  interpolation. 

The  advantages  of  Sauer's  method  are  illustrated  by  a  discussion  of 
its  application  to  the  two  problems  mentioned  above.  The  first  is  most  con¬ 
veniently  referred  to  spherical  polars  and  is  partly  a  finite  difference  method 
while  the  second  is  an  initial  motion  problem  solved  by  expansion  in  series 
and  referred  to  cylindrical  polar  coordinates. 

2.  GENERAL  FORM  OF  SAUER'S  NEAR  CHABACTtBTSTICS  METHOD 

Consider  a  system  of  m  quasilinear  partial  differential  equations  of 


the  first  order 
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k  Sui 
ij  dxi 


2.1 


where 


j  -  1, 


k  -  1,  .. 


1  “  1,  . . . .  n 


and  n  +  1  is  the  number  of  Independent  variables  (Including  the  time  t).  The 
k  lc 

coefficients  a^j,  b^  and  h  are  functions  of  t,  x  and  n. 

We  rewrite  equations  2.1  so  that  only  derivatives  with  respect  to 
t  and  occur  on  the  left  while  all  remaining  derivatives,  together  with 
terms  h  are  put  on  the  right.  We  then  multiply  each  equation  by  a  factor 
CT^  and  add  the  resulting  equations.  We  shall  obtain  the  equation 


du- 


Aih  +  »15il  +  A2^2  +  B2 

dt  Sxj  St  dxj 


2.2 


+  A  +  B 


du 


"St 


-  H 


where 


AJ  “  rk  *kj 


B 
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H  is  a  function  of  (Tm,  t  ,  x^...^  amd  contains  derivatives  with 

respect  to  x2, . xr. 
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We  now  choose  0~.t  ...  G~  so  that  the  directional  derivatives  occurring 
I  in 

on  the  left  of  equation  2.2  are  all  the  same.  The  consnon  direction  is 
defined  by 


±lm\ 

dt  A, 


m 


“  dx  ■  0 

n 


2.3 


Equation  2.3  then  leads  to  m  equations 


®1  '  tA1  “  B2  ‘  tA2 


“  B  -  tA 
m  m 


2.4 


to  determine  ....  ff^.  These  will  give  non-trivial  solutions  only  if 

the  determinant  formed  by  the  coefficients  of  J T  .  <7"  is  zero. 

The  determinantal  condition  gives  r  values  of  tp(p  ■  1,  ...  r).  In  the 
case  of  hyperbolic  equations  2  <  r  <  m  and  the  introduction  of  the  values 
into  equations  2.4  leads  to  m  independent  equations  of  type  2.2.  These 
nay  now  be  written 


di*i 

<-A 


+ 


+  A 


du 


H 

P 


p  *  1 f  ...  r.  k  *  1^  ...  m 


2.4 


Equations  2.4  are  in  every  way  equivalent  to  the  original  equations  2.1  and 
are  transformations  of  these  in  terms  of  near  characteristics.  The  effect 
of  the  use  of  near  characteristics  is  to  reduce  the  number  of  directional 
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derivatives  in  each  of  the  original  equations  by  one,  a  significant 
simplification.  In  the  case  of  one  dimensional  unsteady  motion,  of  course, 
this  transformation  results  in  each  equation  containing  one  directional 
derivative  only. 


3.  UNDERWATER  EXPLOSIONS  IN  AN  OCEAN  WITH  A  VERTICAL  DENSITY  GRADIENT 

This  problem  is  conveniently  referred  to  spherical  polar  coordinates. 
The  Eulerian  equations  of  motion,  for  rotational  flow,  may  be  written 


& 

dt 


+  uS 


+  ± 
p 


w  ou 

7  55 


If 


dw 

dr 


w  dw  uw  _  1  dp 
r  dO  '  r  "  pr 


“tr  +  71  dt  +  oji  dr 
or  pa  ot  pa  or 


a  dw  2ua  wa  .  Q 

- - cot  0 

r  dfi  r  r 


_S_  & 
rpa  dfl 


3.1 


3.2 


3.3 


Here,  t  is  the  time  measured  from  the  Instant  of  initiation,  r  is  distance 
measured  from  the  center  of  the  explosion  and  9  is  the  angle  between  a  vertical 
axis  through  the  center  and  the  radius. 


Multiply  3.1  by  (1^  3.2  by  (T^,  3.3  by  0"^  and  add 

Then 

°i£  ♦  ♦  /St  ♦  <S  +  r>t 


+  <*2  t  *  •*£ 


3.4 
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where  H  is  the  sum  of  the  terms  on  the  right  of  the  factored  equations.  Ihe 
condition  that  the  directional  derivatives  on  the  left  of  3.4  should  all  be 
the  same  is 


0"jU  +  CTjU 


O^a  +  0^u 


<r2u 


Hence, 


(u  -  t)  (T. 


+  a 


-  0 


«  O', 


+  (u  -  t)  0"3  ■  0 


(u  -  t )<r 


-  o 


For  non-trivial  values  of  <7^,  0^,  x  must  satisfy  the  equation 


(u  -  t)  £  (u  -  x)2  -  a2  | 


-  0 


Hence 


3.5 

3.6 

Ihe  equations  of  the  near  characteristics  are  then 


or 


•  u  +  a 


dt 


u  +  a  (Traces  of  Mach  Lines) 


3.7 


8 


*  u  (Projections  of  Particle  Paths)  3.8 

Denote  directional  derivatives  (w.r.t.  time)  in  the  three  characteristic 
conditions  by  suffices  1,  2,  and  3  respectively.  Unis, 


(— ) 
Wl 


&  +  (“  +  a)  ^ 


(-£) 

Mt'2 


£  +  <U  ‘  a) 


t  d*  d  ,  & 

(dt*3  "  at  +  U  ar 


Then  the  near  characteristic  forms  of  equations  3.1,  3.2,  3.3  are  obtained 
by  putting  CTy  ■  1,  0^  ■  +  1,  ■  0;  m  ^3  m  0,  °2  “  1*  We 

obtain 


_  v  _  v  _  a 

vdt'l  pa  at'l  r  r  a®  r  ao 


2a*  .  H*  cot  0  -  -s- 

r  r  rpa  o® 


3.9 


(dU)  .  -±(&) 

'dt'2  pa'dt'2 


fH  +  ? 


2ui 


+  =»  +  “  cot  9  +  & 

r  r  rpa  d9 


3.10 


uw 

r 


9 


,dwv  w  5w 

'dr  3  "r  59 


1  dp 
pr  de 


3.11 


To  start  the  method  we  need  to  know  the  complete  disturbance  field  at 

some  time  t  *  t  .  The  determination  of  these  initial  conditions  will  depend 
o 

on  the  type  of  explosion  considered.  In  general  we  may  suppose  that  they 
are  given  in  the  form 


u  ■  u(r,  9)  etc.  at  t  *  t  3.12 

To  determine  flow  conditions  at  a  later  instant  +  4  t,  at  a  new 
point  Q(rQ,  80)  we  proceed  as  follows:  (1)  Calculate  all  flow  conditions 
at  the  point  P,  which  is  the  projection  of  Q  on  the  t  •  tQ  plane.  In 
particular,  determine  the  characteristic  directions  ^dt^3  at 

P.  (2)  At  Q  draw  lines  with  these  characteristic  directions  to  intersect 
the  plane  t  ■  tQ  at  points  A,  B  and  C  respectively.  Note  that  these  lines 
all  lie  in  the  plane  9  *  0q.  (3)  Write  equations  3.9,  3.10  and  3.11  as 

difference  relations  along  AQ,  BQ,  CQ  respectively.  The  terms  on  the  right 
of  these  equations  are  evaluated  at  A,  B  and  C  respectively.  Solve  these 
equations  to  determine  the  values  of  u,  p,  and  w  at  Q. 

In  the  general  case  these  equations  need  to  be  supplemented  by  the 
equation  of  state  and  the  condition  that  entropy  is  conserved  on  particle 
paths.  The  latter  can  be  written  in  characteristic  form 

.  _  X  M 

W3  r  59  3.13 


which  can  be  solved  as  a  difference  equation  along  CQ.  In  water  the  Tait 
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equation  of  state 


may  be  used.  In  general  B  and.  A  can  be  taken  as  constant  and  equation 
3.13  can  be  dropped. 

This  crude  difference  process  can  be  improved.  Once  first  approximations 
to  the  values  of  flow  quantities  at  Q  have  been  determined  equations  3.9  -  3.12 
can  be  solved  again  along  QA,  QB,  QC  using  mean  differences. 

At  each  stage  of  the  integration  shocks  and  other  boundaries  (e.g. 
contact  discontinuities)  must  be  fitted.  The  procedure  for  doing  this  is 
similar  to  that  used  in  one  dimensional  flow  but  must  be  applied  in  a  series 
of  planes  6  ■  constant. 

Three  sets  of  initial  conditions  can  be  used,  corresponding  to  three 
types  of  explosion.  The  simplest  corresponds  to  a  pressurized  sphere  which 
is  burst  symmetrically  at  time  t  ■  0.  Secondly,  we  may  consider  a  spherical 
charge  of  explosive  initiated  at  its  center;  in  this  case  the  blast  propagation 
is  preceded  by  a  spherically  symmetric  detonation  phase.  The  initial  field 
of  disturbance  in  both  these  cases  can  be  calculated  by  series  expansions 
similar  to  those  used  by  Holt  (ref.  8)  to  analyse  the  initial  behavior  of 
spherically  symmetric  explosions.  The  coefficients  of  these  series  will  now 
depend  on  0,  as  well  as  r  and  t.  For  small  times  after  initiation,  when 
the  vertical  density  gradient  is  small,  it  is  convenient  to  use  Fourier 
series  in  0. 

Thirdly,  nuclear  explosions  should  be  considered.  These  are  determined 
from  perturbation  of  similarity  solutions  and  have  been  treated  extensively  in 
the  book  by  Korobeinikov,  Mel'nikova  and  Ryazanov  (ref.  9). 
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4.  REFLECTION  AND  REFRACTION  OF  A  SPHERICAL  SHOCK  WAVE  AT  THE  OCEAN  SURFACE 


The  problem  of  reflection  of  a  spherical  explosive  wave  at  a  solid  plane 
surface  has  been  solved  by  Vasil' ev  (ref.  10).  Series  expansions  were  used 
in  powers  of  time,  measured  from  the  instant  the  wave  first  hits  the  surface, 
and  radial  distance,  measured  from  the  initial  point  of  contact  of  the  wave. 

By  retaining  terms  in  the  series  up  to  the  fourth  degree  the  shape  of  the 
reflected  shock  wave  and  the  values  of  pressure  and  density  behind  it  can 
be  worked  out  accurately  up  to  the  stage  when  the  incident  shock  makes  an 
angle  of  40°  with  the  plane  surface. 

This  suggests  the  investigation  of  a  more  difficult  problem,  the 
reflection  and  refraction  of  a  spherical  explosive  wave  at  a  free  surface, 
especially  the  surface  of  the  ocean.  If  both  the  ocean  and  air  above  it 
are  treated  as  uniform  media  and  if  a  spherical  explosion  is  detonated  some 
distance  below  the  surface,  the  resulting  wave,  on  reaching  the  surface,  will 
be  partly  transmitted  as  a  shock  into  the  atmosphere  and  partly  reflected 
as  an  expansion  wave  back  into  the  ocean.  In  any  meridian  plane  this  expansion 
wave,  at  a  given  time,  is  centered  on  the  point  in  the  surface  then  reached 
by  the  incident  wave.  The  ocean  surface  is  disturbed  by  the  refraction  process 
and  it  is  necessary  to  determine  the  shape  of  the  transmitted  wave,  the  dis¬ 
tortion  of  the  ocean  surface,  and  the  field  of  interaction  between  the 
reflected  expansion  wave  and  the  incident  explosive  wave  (see  Figure  2). 

To  treat  this  initial  motion  the  governing  equations  are  written  in 
cylindrical  polars. 

Momentum 


St  +  "Sr  p  Sr 


4.1 


12 


& 

at 


+ 


(to  _  1  dp 

dr  p  dz 


4.2 


Continuity 


& 

St 


+ 


P* 


& 

dz 


4.3 


Sat ropy  conservation 


& 

St 


+ 


4.4 


State 


B(S) 


£  -  A(S)  | 


4.5 


In  water  B  and  A  are  taken  to  be  constants;  y  ■  7 
In  air  B  -  e^^v,  A  -  0,  Y  -  1.4 

The  z  axis  Is  taken  on  the  vertical  line  through  the  origin  of  the 
explosion.  The  radial  coordinate  r  Is  measured  from  the  point  on  the  ocean 
surface  where  the  spherical  shock  is  first  in  contact. 

To  determine  the  field  of  disturbance  near  the  point  ( r ,  0)  it  is 
convenient  to  introduce  a  reduced  time  t.  This  is  measured  from  the  instant 
when  the  incident  wave  reaches  the  point  (r,  0)  rather  than  the  instant  of 
initial  contact. 

If  U(&)  is  the  value  of  the  incident  shock  wave  velocity  (K  is  measured 
from  the  center  of  the  explosion)  the  point  of  intersection  of  the  incident 
shock  with  the  free  surface  is  moving  with  velocity 

°1  "  “3t  "  h  [  °<h>  +  2  <hD'<h>  -  0(h))  ...j  4.6 
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where  h  is  the  depth  of  the  center  of  the  explosion 
The  reduced  time  t  is  therefore , 


T 


t  - 


t 


1 


or 


+  \  “5 -  (hU'(h)  -  0(h))  ... 

2  MT(h) 


4.7 


In  equations  4.1  -  4.4  change  independent  variables  from  r,  z,  t  to  r,  z,  t. 
Then 


where 


_S 

Sr 


kr 


-3 


St 


(neglecting  higher  order  terms  in  r) 


h  U^h) 


" 


St 


. 

S* 


.a 


S* 


4.8 


4.9 


4.10 


4.11 


The  transformed  equations  are  to  be  expanded  in  series  of  powers  of 

r.  To  determine  the  lowest  order  terms  in  these  expansions  we  may  neglect 
2 

terms  of  order  r  .  Now  u  is  of  order  r  near  r  ■  0  and  we  may  therefore 
neglect  kru  (and  kr  ^“). 
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The  equations  become 


Su  .  chi  ,  kr  op 

+  +  rsT 


-  u* 


Su 

Sr 


p  tr 


&  +  w  &  +  I  &  .  .  „^Si 
dr  dz  p  Sz  uSr 


Sf 


,  Sp  2  dw 

+  w  +  pa  — 


dz 


9z 


2  ,au  ,  uN 
-  pa  +  ~) 


as  .  as 

f  +  wa? 


-  u 


as 

ar 


He  now  write  equations  4.12  -  4.14  in  near  characteristic  form.  The 
linear  combination  of  equations  4.12  -  4.14  is 


ci$  ♦  VS  +  <*-3  +  ht>£  -  %  *  V>l 
+  r2  £  *  <  v +  v2>  S 
■  -  <V  +  V'**  -  ^  t  -  V  £  -  V*2  ? 


The  near  characteristic  conditions  are  determined  by  the  conditions 


<r2 

££  .  ~+  r3w 

<rx  "  <r3  +  o^kr 

p 


°2W  +  ^3P* 


1 


4.12 


4.13 


4.14 


4.15 


general 


4.16 


15 


Hence 

(w  -  X)  CT  ■  0 

-  —  x  cn  +  —  +  (w-x.)cr  -  o 

pip  3 

(w  -  X)  <7"2  +  pa2  0~3  *  0 


For  non-trivial  <7^,  0^ 

(w  -  X)  0  0 

-  1  (w  -  X> 

P  P 

0  (w  -  X)  pa^ 


-  0 


Hence 


X  ■  w  (particle  path  direction) 

or  X  ■  w  +  a  (ray  directions)  4.17 

With  w  take  0^  ■  1,  £7^  -  krw,  0^-0 

Then  equation  4.16  gives 


(— ) 
Mt;3 


k£ 

P 


<i£) 

'dt'3 


k™  <J?>3 


-■fc 


.  1  £& 

p  Sr 


.  Sw 

krwu  Sr 


4.18 
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where 


'dt'3 


w 


_a 

Sz 


With  \  m  w  +  a  take  0^  ■  0,  O'  -  1,  O'  *  +  1/ (pa) 
Equation  4.16  then  reduces  to 


+  ~  (■ 
-  pa 


dv 

dt/l,2P 


+  dr) 


.  dx  _  —  Su  Sw  —  au 

■Tr),  -  +  a?*  -  u?=  +  — 


dt  1,2 


Sr  “Sr 


where 


'dt^l,2 


(w  +  a) 


_S 

Sz 


Equation  4.15  can  be  written 


(*>) 

Mt'3 


4.19 


4.20 


4.21 


4.22 


Equations  4.18,  4.20,  4.22  are  in  the  required  near  characteristic 
form.  To  determine  the  initial  flow  pattern  shortly  after  the  incident 
wave  first  hits  the  ocean  surface  we  expand  the  dependent  variables  in  powers 
of  r,  for  example. 


P  "  PQ(t»  *>  +  r  Pi(t,  *> 


and  retain  terms  up  to  order  r.  We  then  obtain  equations  for  the  coefficients 
written  in  characteristic  form  and  these  are  solved  in  series  of  the  type 

po  -  Poo^  +  pol(f> 
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starting  with  the  zero  order  coefficients.  These  expansions  are  very  similar 
to  those  carried  out  for  the  spherical  explosion  by  Holt  (ref.  8). 

The  detailed  calculations  for  both  these  problems  are  now  being 
carried  out  and  results  will  be  published  in  later  reports. 
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Figure  1. 


t 


A  C  B 


Difference  Scheme  in  Hear  Characteristics  Method 


Figure  2.  Reflection  and  refraction  of  a  spherical  shock  wave. 
ACA'  Transmitted  shock. 

ADA1  Disturbed  ocean  surface. 

AEA'  -  AFA'  Reflected  expansion  wave. 

AB,  A'B'  Incident  shock. 


